Analytical and numerical study of trapped strongly 
correlated bosons in two- and three-dimensional lattices 
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We study the ground-state properties of trapped inhomogeneous systems of hardcore bosons in 
two- and three-dimensional lattices. We obtain results both numerically, using quantum Monte Carlo 
techniques, and via several analytical approximation schemes, such as the Gutzwiller-mean-field 
approach, a cluster-mean-field method, and a spin-wave analysis which takes quantum fluctuations 
into account. We first study the homogeneous case, for which simple analytical expressions are 
obtained for all observables of interest, and compare the results with the numerical ones. We obtain 
the equation of state of the system along with other thermodynamic properties such as the free 
energy, kinetic energy, superfiuid density, condensate density and compressibility. In the presence 
of a trap, there is in general a spatial coexistence of superfiuid and insulating domains. We show 
that the spin-wave-based method reproduces the quantum Monte-Carlo results for global as well as 
for local quantities with a high degree of accuracy. We also discuss the validity of the local density 
approximation. Our analysis can be used to describe bosons in optical lattices where the onsite 
interaction U is much larger than the hopping amplitude t. 

PACS numbers: 03.75.Hh,03.75.Lm,67.85.-d,02.70.Ss 
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I. INTRODUCTION 

The successful realization of the superfluid-to-Mott- 
insulator transition in ultracold bosonic gases trapped 
in optical lattices in one [1| , two @, Q , and three [J] di- 
mensions has opened the way to a myriad of experimen- 
tal and theoretical studies of strongly correlated lattice 
systems @ . One of the main ideas driving these stud- 
ies is that ultracold atoms in optical lattices can be used 
as analogue simulators of Hamiltonians of the Fermi- and 
Bosc-Hubbard type. Intensive efforts arc currently under 
way to validate this approach by comparing experimen- 
tal and theoretical results for systems such as the Bose- 
Hubbard model which is amenable to both treatments 

In the Bose-Hubbard model, the basic Hamiltonian 
consists of a hopping term of amplitude t and an onsite 
two-body repulsion term with amplitude U . The phase 
diagram of this homogeneous model is known to consist 
of two phases: (i) a superfiuid phase that is present for 
all incommensurate fillings and arbitrary values of the 
ratio U /t, and for commensurate fillings below some crit- 
ical value (U/t)c, which depends on the dimensionality of 
the system and on the (integer) filling, and (ii) a Mott 
insulator which is p resent for commensurate fillings for 

u/t > {u/t)e Um- 

In experiments with ultracold gases, a confining po- 
tential is always present. This confining potential is to 
a good approximation harmonic, and generates an inho- 
mogeneous density profile in which superfiuid and Mott- 
insulating phases coexist in spatially-separated domains 
|13l - [l8j . In such systems, the appearance of Mott insu- 
lating domains in different regions of the trap depends 
not only on the total filling and the ratio U/t (as in the 
homogeneous case) but also on the curvature of the har- 



monic confining potential. 

So far, an accurate characterization of the local prop- 
erties of the system in the trap has only been achieved 
by means of quantum Monte Carlo (QMC) simulations 
[lj,[l5|,[l3,[l3: and density-matrix renormalization group 
(DMRG) techniques in one dimension [la|. Our goal in 
this paper is to introduce an analytical approach that 
provides an accurate prediction of the behavior of the sys- 
tem in two and three dimensions. In this study we shall 
consider the strongly correlated limit where the ratio U/t 
is very large and the local occupation of the lattice sites is 
lower than or equal to one. In that case, the system can 
be described, to a good approximation, by impenetrable 
(hardcore) bosons. As we will show, this in turn enables 
the implementation of an analytical treatment which goes 
beyond simple mean-field calculations, where within this 
approach, quantum fluctuations are taken into account 
by the addition of spin-wave corrections. 

In what follows, we shall examine the extent to which 
the results of this method arc valid by comparing them 
against quantum Monte Carlo simulations. As we will 
show, spin- wave-corrected results provide a big improve- 
ment over the usual Gutzwiller-mean-field approach and 
its cluster-mean-field extension. In most cases, for the 
system sizes considered here, they are also more accu- 
rate than the QMC-based local density approximation 
(LDA). In some regimes, and for some local observables, 
the latter is found to be less accurate than the simple 
Gutzwiller-mean-field approach. 

The paper is organized as follows. In Sec. |TT] we 
review the model at hand and discuss its properties. 
In Sec. mil we study homogeneous systems in two and 
three dimensions. We obtain the equation of state along 
with basic thermodynamic properties and examine the 
extent to which several analytical approaches compare 



against exact numerical QMC simulations. In Sec. IIVI 
we analyze the global and local properties of hardcore 
bosons trapped in a confining harmonic potential in a 
two-dimensional setup. We also examine the validity of 
the local-density approximation for finite systems and 
compare the various analytical approximations against 
exact QMC simulations. In Sec.|Vl we analyze the behav- 
ior of three-dimensional harmonically trapped hardcore 
bosons using the various approximation schemes. Finally, 
in Sec. I VII we conclude with a discussion and summary 
of our results. 



II. MODEL 

The Hamiltonian for hardcore bosons confined in a 
harmonic potential in a d-dimensional hypercubic lattice 
with N = L'^ sites, can be written as: 



H 






+ a] a, 



Ai 



H "' 



vY^r^n,. (I) 



Here, (ij) denotes nearest neighbors and hi {a\) destroys 
(creates) a hardcore boson on site i located at a distance 
ri = \ri\ from the center of the trap, where the coordi- 
nates Ti = {xii, . . . ,Xdi) are given here in units of the 
lattice spacing a, which we set to unity. 

The operator fii = a\ai is the local density operator, /i 
is the global chemical potential, and V is the harmonic 
potential strength. The hopping parameter t (which we 
shall fix at i = 1) sets the energy scale. The hardcore bo- 
son creation and annihilation operators satisfy the con- 
straints a\ = a| = 0, which prohibit double or higher 
occupancy of lattice sites, as dictated by the C/ — > oo limit 
of the Bose-Hubbard model. For any two different sites 
i ^ 2i the creation and annihilation operators obey the 



usual bosonic relations 



0.1, a j\ 



\a\,a\\ 



\ai,a\ 



0. 



To gain a general understanding of the zero- 
temperature phases of this model, let us first analyze 
the atomic (t — 0) limit. In this limit, there is no ki- 
netic (hopping) term, and the boson number operators 
hi commute individually with the Hamiltonian, so every 
lattice site is occupied by a fixed number of bosons. In 
this case, the Hamiltonian is diagonal in the Fock-states 
basis, and the ground-state wave-function has the form 
|GS) = Jlf" iV'i) and is solved for each site separately. 



givmg 



\i'^) 



|0) if /i, < 
|1) if Aiz > 



(2) 



where we have denoted yi,i = [p. — Vrf) as the local chem- 
ical potential, and |0) and |1) denote vacant and occupied 
sites, respectively. In the special V = ^ case, where no 
trap is present, the model is translationally invariant, and 
the ground-state boson occupancy is the same through- 
out the lattice: For /^ < the minimal energy configu- 
ration is simply the vacuum, i.e., the completely-empty 



lattice, and for /x > the minimal energy configuration 
is the completely-filled lattice. The ground-state energy 
of these phases is degenerate at /i = 0. 

In the t ^ case, the model has in general no analytic 
solution [13 • As it turns out, however, the phase bound- 
aries separating the insulating phases from the superfluid 
one can easily be obtained analytically in the 1/ = case 
even for nonzero t. To see this, we use the fact that our 
Hamiltonian commutes with the total-number-of-bosons 
operator iV;, = X^i^*- This simply means that for any 
given /i and t, the ground-state wave-function is a lin- 
ear combination of product states each having the same 
number of occupied (vacant) sites. In the completely- 
filled phase (which we denote by F) the wave-function is 
simply 



iGS)F=nii) 



(3) 



with energy ep = —fiN. In the infinitesimally thin layer 
outside of the F phase, the state of the system is char- 
acterized by exactly one vacant site, and thus its wave- 
function is the sum 



|GS)f+ =A^"^/2^a,;|GS)j. 



(4) 



with energy ep+ — —2dt — ^{N — 1). The phase bound- 
ary separating the superfluid phase and the (insulating) 
completely filled phase is the curve along which the F 
state, Eq. ([3]), is no longer energetically favorable. This 
happens when its energy becomes degenerate with the 
energy of the defect state, Eq. ^. Matching the two, we 
obtain the phase boundary: 



2dt 



1, 



(5) 



By the same token, the boundary between the superfluid 
phase and the completely empty insulating phase (which 
we denote by E) is given by fi/{2dt) = — 1. This result 
can also be obtained by repeating the above exercise with 
the substitution |0) O |1). Between the two insulating 
phases, in the region \^/2dt\ < 1, the system is super- 
fluid. 

When a trap is introduced, the system will no longer 
be entirely superfluid or entirely insulating. Depending 
on the trap curvature and the total chemical potential, 
the system may be in a coexistent state of a superfluid in 
one part of the lattice and an insulator in another. Before 
discussing why this is so, let us flrst note that unlike the 
homogeneous case which is scale- invariant, in the pres- 
ence of a trap the model has a length scale determined by 
the ratio of the trap curvature to the hopping amplitude, 
given by ^ = {V/t)~^''^. As a result, our Hamiltonian 
can be rescaled by introducing the dimcnsionless scaled 
length fi = Ti/S, through which the Hamiltonian can be 
rewritten in a way that eliminates the amplitude of the 



harmonic trap V: 



H = -t 






t 



-f2 



(6) 



As a result of this rescahng, it is clear that systems with 
the same global chemical potential exhibit the same qual- 
itative behavior for all positive trap curvatures. In addi- 
tion, away from the center of the trap, the local onsite 
potential fii always becomes very large and negative and 
hence outside of some radius fout, the lattice must be 
empty {{fii) = 0). On the other hand, in the center of 
the trap, large enough local chemical potentials will pro- 
duce an insulating (circular) domain with (hi) = 1 inside 
some critical radius r < fjn. In the intermediate region, 
Hn l£ r < fout, the system is in the so-called superfluid 
regime. The radial density profiles of the two possible 
scenarios, with and without the insulator in the center of 
the trap, are sketched in Fig. [T] 
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FIG. 1: Typical density profiles of harmonically trapped hard- 
core bosons as a function of the distance from the center of 
the trap. For large enough values of chemical potential (solid 
line), a Mott insulator with (hi) = 1 forms in the center of the 



trap (the region f < rlJ) 
is superfluid up until f — 



Outside of this radius the system 
outside of which the lattice is 



d^) 



empty. If the chemical potential is small (dashed line) a Mott 
insulator will not form in the center of the trap, and the only 
superfluid- insulator transition will take place at f = fout- 



III. HOMOGENEOUS (V = 0) CASE IN TWO 
AND THREE DIMENSIONS 



Having discussed the general properties of the model 
in the previous section, we now turn to analyze it in a 
quantitative manner. We shall try to keep the discussion 
as general as possible so that it applies to arbitrary di- 
mensions and confining potentials. Only after the general 
formulation has been introduced, we will obtain analyt- 
ical expressions for homogeneous (i.e., V = 0) systems 
in two and three dimensions. (We note that some of the 
observables of the homogeneous two-dimensional case we 
analyze here have been studied in Ref. [201.) This general 
formulation will be useful later on when we address the 
problem of harmonically-confined bosons. 



Our main objective in this section is to obtain the equa- 
tion of state of the model, i.e., the average density as a 
function of the chemical potential, along with the basic 
thermodynamic properties of the system. Since in dimen- 
sions higher than one the inodel has no analytic solution, 
exact results for finite systems (up to the relevant statis- 
tical errors) arc obtained numerically. Here we use the 
stochastic series expansion (SSE) algorithm [2l|, ^^ and 
perform simulations over a wide range of chemical po- 
tentials. Zero-temperature properties of the system are 
obtained by choosing large inverse-temperatures (3 = 1/T 
(in our units, ks = 1), where we have found it sufficient 
to have /3 > L, as the effects of increasing /3 beyond this 
value are indiscernible. In two dimensions, we simulate 
systems with 48 x 48 sites and 60 x 60 sites, with an 
inverse-temperature of /3 = 60, and in three dimensions, 
simulations arc performed on a 16 x 16 x 16 lattice and 
f3 = 20. After obtaining the QMC results, we examine 
the validity of several approximation schemes by compar- 
ing them against the results of the SSE simulations. 



A. Quantum Monte Carlo results 

Wc first study the equation of state, which is very rele- 
vant for experiments as it provides the dependence of the 
number of bosons in the system (equivalently, the average 
density) on the global chemical potential. It is depicted 
in Figs, m^ a) and[31Ja) for two- and three-dimensional lat- 
tices, respectively. The corresponding compressibilities 
defined by k = dp/dfj, are also given, in Figs. [2Jb) and 
^h). These provide the response of the average density 
to a change in the chemical potential. We have com- 
puted them in two independent manners: (i) as the nu- 
merical derivative of the average density with respect to 
the chemical potential and (ii) within the SSE algorithm 
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FIG. 2: (Color online) (a) Equation of state and (b) com- 
pressibility K = dp/djj, for hardcore bosons in a homogeneous 
two-dimensional lattice (60x60 sites). The SSE results are de- 
noted by circles, whereas the dotted, dashed, and solid lines, 
are the results of the Gutzwiller- mean-field, the cluster-mean- 
field, and the spin- wave-corrected Gutzwiller-mean-field ap- 
proximations, respectively. 



using the formula: 
_dp _ 1 



iVfcC 



-I3H 



^j^[m)-m') ,(7) 



where Z is the partition function and we have taken ad- 
vantage of the fact that [Ni,,H] = 0. As expected, the 
results of both methods were found to coincide. In the 
figures, the SSE results are denoted by circles, whereas 
the various lines indicate the predictions obtained by sev- 
eral approximation schemes that will be introduced later. 
As the figures indicate, at fi/{2dt) = the average den- 
sity is 1/2. This is because of the particle-hole symmetry 
of the model. Upon increasing the chemical potential, 
the boson density increases too until it reaches p = 1 at 
^/(2dt) = 1, as predicted by the analysis given in the 
previous section. At p = 1 the system becomes insulat- 
ing and the compressibilities drop sharply to zero. Note 
that due to the particle-hole symmetry, the results for 
negative fi can be immediately read off the figures by the 
transformations p, ^f —p, and p — >■ 1 — p. 



0.05 




FIG. 3: (Color online) (a) Equation of state and (b) com- 
pressibility K = dp/ dp for hardcore bosons in a homogeneous 
three-dimensional lattice (16 x 16 x 16 sites). The SSE re- 
sults are denoted by circles, whereas the dotted, dashed, and 
solid lines are the results of the Gutzwiller-mean-field, the 
cluster-mean-field, and the spin-wave-corrected Gutzwiller- 
mean-field approximations, respectively. 

Figures 2] and [5] show other observables of interest in 
two and three dimensions, respectively. Those observ- 
ables provide further understanding of the properties of 
the model: Figures Slja) and ^a) show the free energy 
(per site), which is the quantity that is being minimized 
throughout the simulations. Also shown [Figs. Hl^b) and 
[SJb)] is the kinetic energy (per site), which can in prin- 
ciple also be measured in ultracold gases experiments by 
time of flight expansion after simultaneously switching 
off the trapping potential and the lattice. 

The superfluid densities [Figs. |4|c) and ^c)] and the 
condensate densities [Figs.Ul^d) andjSJd)] further indicate 
that for p < 1; when the systems are compressible, they 
are also superfluid and exhibit Bose-Einstein condensa- 
tion. In the presence of the strong correlations generated 
by the infinite onsite repulsion, one can see that the su- 
perfluid density is always greater than the condensate 
density. This is similar to what happens in liquid helium 



where strong correlations deplete the condensate density 
to a very small value while the superfluid density remains 
very large. Interestingly, as the dimensionality of the sys- 
tem increases one can see that the difference between the 
superfluid density and the condensate density decreases. 
This point will be touched upon later, when we discuss 
the Gutzwiller-mean-field method which corresponds to 
the exact solution in infinite dimensions. Also evident 
from the figures is the fact that the superfluid density 
and the condensate density become equal in the low den- 
sity limit. 



B. Approximation schemes 



Having analyzed the homogeneous case via the SSE 
algorithm, we now proceed to analyze the model with 
a number of approximation schemes. We start this in- 
vestigation with the Gutzwiller mean-fleld approach. We 
note here that this approach is a particular case of a more 
general treatment initially developed for soft-core bosons 

S0- 



0.5 



P 

0.75 



10.5 



P 
0.75 




0.26 



.&■ 



6 0.13 





-0.5 



0.24 



'£- 



-0.12 



0.5 



0.75 
P 



1 0.5 



0.75 
P 



FIG. 4: (Color online) Thermodynamic properties of hard- 
core bosons in a homogeneous two-dimensional lattice: (a) 
free energy, (b) kinetic energy, (c) superfluid density, and 
(d) condensate density as a function of the average den- 
sity. The circles indicate the SSE results (60 x 60 sites), 
whereas the dotted, dashed, and solid lines indicate the 
Gutzwiller-mean-field, the cluster-mean-field, and the spin- 
wave-corrected Gutzwiller-mean-field approximations, respec- 
tively. 
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FIG. 5: (Color online) Thermodynamic properties of hard- 
core bosons in a homogeneous three-dimensional lattice: (a) 
free energy, (b) kinetic energy, (c) superfluid density, and 
(d) condensate density as a function of the average den- 
sity. The circles indicate the SSE results (16 x 16 x 16 
sites), whereas the dotted, dashed, and solid lines indicate the 
Gutzwiller-mean-field, the cluster-mean-field, and the spin- 
wave-corrected Gutzwiller-mean-field approximations, respec- 
tively. 



1. Gutzwiller mean-field 

Generalizing the discussion in Refs. [23 and [2J1 to 
inhomogeneous systems, we start our mean-field calcula- 
tion with the following product state as our ansatz for 
the ground state wave-function: 



|GS) 



MF 



n 



sin^lO), 



cos^e^'^^ll), 



(8) 



where the unknowns {9j , Lpj ) are to be determined via the 
minimization of the grand-canonical potential (per site) 
which evaluates to 



r^MF = 



1 



MF 



(GS|iJ|GS)MF 



(9) 



— — ^ sin 6*4 sin 9j cos (</>,; - 0j) 



2N 



(u> 



" T^^t^^O^+^^^^i) ■ 



In the homogeneous case, the wave- functions of each of 
the lattice sites are identical. However, in the following, 
we consider general local potentials, as later we shall ap- 
ply this method to the case of the harmonic potential. 
For the azimuthal angles ipj^ it is sufficient to require 



that all sites share the same value (pj = $, where $ can 
be chosen arbitrarily. The conditions for the polar angles 
are obtained by differentiating JImf with respect to 9j, 



giving: 



fii tan t 



ts. 



(10) 



J — -j"i,n(i) wnere Oj,„(j) is unity if i and j 



Here, Si = J2j sm(^j^i.n{j) where S, 

are neighbors and is zero otherwise. This set of equations 

is to be solved numerically in the general case. 

In the limit where the lattice spacing is very small com- 
pared to the average interparticle distance, i.e., in the low 
density limit, the ^^'s may be assumed to vary smoothly 
over the lattice and the above set of equations may be 
approximated by its continuous version, which turns out 
to be the ordinary differential equation 



fi{r) tan 0{r) = tV^ [sin 6l(r)] + 2dt sm9{r) 



(11) 



which also has to solved numerically in the general case. 
In the homogeneous case, where V = and the sys- 
tem is translationally invariant, Eq. (|10p can be solved 
exactly, as we can write dj = 6 from which it follows that 



cos 9 = Min 



l,Max 



-1 ^ 
' 2dt. 



(12) 



Now that the ground state wave-function has been ob- 
tained, all physical properties of interest can easily be 
calculated by computing the expectation value of the ap- 
propriate operator. The average density of particles, for 
example, becomes 



Pmf 



— 2_^MF(GS|ajai|GS)MF 

i 

^E(i 



cos9t) , 



(13) 



where in the homogeneous case, this expression simplifies 
to 



1 , 
Pmf = - (1 + cos( 



(14) 



The density of bosons is plotted in Figs. [2] and [3] for two 
and three dimensions, respectively, along with the cor- 
responding compressibilities, as a function of the chem- 
ical potential. As can be seen in Figs, ^a.) and[3ja), 
when compared against the quantum Monte Carlo re- 
sults, the Gutzwiller mean-field approach provides a good 
approximation for the density close to half filling and 
again at ^/{2dt) — 1 where the lattice becomes com- 
pletely filled, but deviates from the QMC results in be- 
tween. The behavior of the compressibility, as predicted 
by the Gutzwiller mean-field approach, is on the other 
hand qualitatively incorrect for all densities in the super- 
fluid regime. 

Next, the free energy, Eq. ^, becomes in the homo- 
geneous case 



n 



MF 



N 



MF 



(GS|F|GS) 



MF 



dt , ^ I , 



(15) 



and the density of bosons in the zero-momentum mode 
(the condensate density) turns out to be 

PO.MF = "77MF(GS|a^^^aj^^p|GS)MF 



mT.^' 



4iV2 



sm Ui sm t 



-sin^e, (16) 



where a, (^fc=o) creates (destroys) a particle with mo- 
mentum k. The superfluid density, Ps.mf, requires a spe- 
cial treatment of the boundary conditions. As is well 
known [231 ' the superfluid density is related to the "spin 
stiffness" of the system: To accomplish this, one needs 
to compare O (the free energy) of the system under pe- 
riodic conditions with the free energy under a "twist" 
in the boundary conditions along one of the linear di- 
rections (say, the x direction). In the homogeneous case 
we consider here, the azimuthal angles ipj are all iden- 
tical. To implement a twist, we take this angle to be 
site-dependent and with a constant gradient such that 
the total twist across the system in the x direction is tt, 
namely Sip ~ fj+x — ^j ~ tt/ L. Within the mean-field 
treatment, one can show that addition of this gradient 
is equivalent to substituting t -^ t/d[{d~ \) + cos&ip]. 
Now, the square of the gradient twist is related to the 
superfluid density via the relation Jltwisted — O = tpsSip'^ 
which in turn yields the simple expression 



Ps,MF 



1 dn 

2d~dt 



(17) 



At the mean- field level, this expression evaluates to that 
of the condensate density. Since mean-field theory be- 
comes exact in infinite dimensions, the differences be- 
tween both observables should decrease as the dimension 
of the systems increases. This explains why the differ- 
ences between the superfluid density and the condensate 
density (as they were computed with the QMC simula- 
tions in Sec. lIII.^ were smaller in three dimensions than 
they were in two dimensions. 

The quantities discussed in Eqs. ([T5|) - pT| are plotted 
in Figs. 2] and [5] for two and three dimensions, respec- 
tively. As one can immediately see, the Gutzwillcr-mean- 
field predictions agree rather poorly with the QMC data 
with respect to all quantities but the free energy, particu- 
larly close to half-filling. Hence, the mean-field approach 
described in this section does not provide an accurate pic- 
ture for most observables of interest close to half filling. 
Nonetheless, it is clear that the agreement with the QMC 
data is better in three dimensions: As already noted, 
the mean-field approximation becomes exact in the limit 
where d — ;> oo. 



2. Cluster mean-field 

The cluster-mean-field approach is an improvement 
over the simple Gutzwiller-mean-field approach, which is 



also applicable to harmonically trapped hardcore bosons. 
It was introduced in Refs. [23| and Q to improve the 
Gutzwiller-mean-field prediction for the phase diagram 
of hardcore bosons when a period-two superlattice was 
added to the homogeneous model. Within this approach, 
one starts with a variational ansatz which is also a prod- 
uct state, but not of single-site wave-functions. The new 
ansatz is a product of wave-functions each describing the 
state of a 'block' of 2'^ sites. In two dimensions, a block 
consists of 2 X 2 square cells, so the ground state wave- 
function has the form: 



|GS) 



CMP — 




-^ijkl 



\ljkl) 



(18) 



Analogously, in three dimensions, the basic block 
is a 2 X 2 X 2 cubic cell ^^. As with the 
Gutzwiller-mean-field case, we minimize the free energy 



CMF 



efficients 



N cMF(GS|i/|GS)cMF with respect to the co- 
of the wave-function. Obtaining the var- 



'ijkl 



ious observables in terms of the wave-function given in 
Eq. (jlSp is straightforward, and is performed in much 
the same way as with the usual mean-field approach dis- 
cussed in Sec lIIIffT] 

The two- and three-dimensional equations of state, as 
determined by the cluster mean-field approximation, are 
depicted by the dashed lines in Figs.[2]and[3l respectively. 
The predictions for the other observables are presented 
in Figs, m and [5] As the figures indicate, the cluster- 
mean-field predictions are an improvement over the plain 
Gutzwiller-mean-field approach but are still far from ac- 
curately reproducing the results of the QMC simulations. 



3. Addition of spin-wave corrections 

In what follows, we show that the results of the mean- 
field approaches discussed above can be significantly im- 
proved by the addition of spin-wave corrections which 
take into account quantum fluctuations J26l - |29{ . This 
method was also discussed in the context of homogeneous 
two-dimensional systems in Ref. |2C| . In a later study 
[23j . it was shown that in the presence of a superlat- 
tice, spin wave corrections do not modify the predictions 
of the simple Gutzwiller-mean-field theory for the phase 
diagram. Here, we generalize this method to arbitrary 
dimensions and general inhomogeneous systems. 

The spin-wave analysis begins with introducing a set 
of local rotations to the boson creation and annihilation 
operators. This is accomplished by switching to new field 
operators via the transformation: 
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The new annihilation and creation operators bj and b- 
describe low-energy fluctuations about the mean-field 
ground state - these are spin waves. They too obey 
hardcore bosons commutation relations. Substituting 
these expressions into our Hamiltonian, ignoring cubic 
and quartic terms in these bosonic operators (thus as- 
suming a dilute gas of spin waves) , the new Hamiltonian 
becomes a sum of a constant term Hq , a linear term Hi 
and a quadratic term H2, where 
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and the matrix coefficients in Jfo are: 



Aij = {^fJ-i COS 9i + t sin 9iSi) 5, 
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(21) 



The spin- waves analysis continues with the determina- 
tion of the 9i^s via the requirement that the linear term 
Hi, Eq. (pn|) . vanishes. Equating each term in the sum 
to zero, one ends up with the set of equations: 



m tan 9i = -tSi , 



(22) 



which is analogous to the set of equations obtained previ- 
ously when the Gutzwiller-mean-field approach was dis- 
cussed, and which similarly has to be solved numerically 
in the general case. In the homogeneous case where 
Hi = fi, we have 9i = 9 and therefore also Si = 2ds[ii9. 
This leads to the solution: 
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Substituting the solution into the zeroth order Hamilto- 
nian Hq, Eq. pO|) . the zeroth order free energy yields, 
as it should, the free energy of the Gutzwiller mean-field 
approach: 
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In the homogeneous case, the above expression evaluates 
to: 
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which is precisely the mean- field result, Eq. ((T5|) . 

In order to obtain the spin- wave corrections, one must 
diagonalize the remaining quadratic term H2 in Eq. (j20p , 
for which one assumes that the spin waves are very dilute 
and hence that the hardcore constraint can be ignored. 
The diagonalization process is a lengthy but straight- 
forward process. Here, we shall only review the basic 
steps. The interested reader is referred to Appendix 
B of Ref. [231 fo'^ ^ more detailed account. The diag- 
onalization process consists of the following steps: (i) 
Diagonalize the matrix {A + B){A — B) and obtain 
its (non-negative) eigenvalues A^ along with the unnor- 
malized eigenstates (/>j, = {4>ki, ■ ■ ■ , 't'ki, ■ ■ ■ , 0feAf} where 
k ^ I . . .N (here, A and B are matrices whose elements 
are Aij and Bij, respectively), (ii) Evaluate the set of 
vectors xpf, = {ipkij ■ ■ ■ ,tpkh ■ ■ ■ , '4'kN} through the rela- 
tion {A — B)4)^. = ^kfpk ■ (iii) Normalize (f}/. and -0^. 
according to J^k '^k ' '4'k = 1- (i'^) Obtain the basic two- 
particle operators through the inverses of cpf. and ^pl.: 

(blbl) = I E i^M ^u) (<^™1 + ^„}) . (26) 



(v) Express any desired physical quantity in terms of the 
new field operators using the transformation, Eq. (fT9l) . 
neglecting cubic and quartic terms, and use the expres- 
sions for the two-particle operators, Eqs. ([25)1 . to evaluate 
the resulting expression. 

As an example, consider the local density of particles. 
Within the spin- wave analysis, it evaluates to 

P^i = (ala,) = ^ (1 - cos 6*0 + cos9^{blbi) , (27) 
following which the average density becomes: 

-"sw = ^E <«!"') (28) 
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The potential energy (per site) is then: 
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The off-diagonal terms are more cumbersome to calculate 
but can still be obtained in a straightforward manner. 
They evaluate to: 
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Once these are calculated, one can obtain the kinetic en- 
ergy (per site): 
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and also the momentum distribution function defined by: 
n(fc) = ^E^"»°"^''''"'"'"'"^- (32) 



In the homogeneous case, the above equations simplify 
substantially, and one can obtain analytic expressions for 
the various physical observables 123| . The spin- wave cor- 
rected density of particles is given by: 
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where we have defined 
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The spin- wave corrected density is plotted in Figs. [2] 
and [3] for the two- and three-dimensional cases, respec- 
tively. As is evident from the two figures, the spin- wave 
corrected densities provide a decisive improvement over 
the plain Gutzwiller-mean-field results and they are also 
closer to the exact results than the cluster mean-field cal- 
culations for the most part. Interestingly, we find that 
for the compressibility, the results of the cluster mean- 
field approach are slightly better than those provided by 
the spin-wave corrections in three dimensions. 

The spin-wave corrections also yield the following ex- 
pression for the free energy: 
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where the superfluid density is obtained by evaluating 
the expression Ps,sw = — (2d)~^ 9risw/i9i. Finally, the 
density of bosons in the zero-momentum mode is found 
to be: 



Po,sw = PO,MF 
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The above quantities are plotted in Figs. |4] (two dimen- 
sions) and [S] (three dimensions). As is clear from the 
figures, not only are the spin-wave corrected predictions 
a major improvement over the mean-field approaches for 
the thermodynamic quantities, but they are essentially 
on top of the QMC results for the kinetic energy and the 
superfluid density, and very close to the QMC results for 
the free energy and the condensate density. 



IV. TRAPPED HARDCORE BOSON IN TWO 
DIMENSIONS 

Having discussed the properties of the homoge- 
neous case, we are now ready to analyze the more 
experimentally-relevant case of lattice bosons in the pres- 
ence of a harmonic trap. We shall limit the discussion to 
the two-dimensional case in this section and address the 
three-dimensional case in the next section. 

As discussed in Sec. |lTl in the presence of a harmonic 
potential, the trapping amplitude V can be scaled out 
by using the dimensionless parameter ^ = {V/t)~^''^. 
In order to verify that this is indeed the case, we have 
performed simulations over a range of chemical poten- 
tials for two different trapping amplitudes, V = 0.01 and 
V = 0.02. We have verified that once properly normal- 
ized, results for local quantities for the two trapping po- 
tentials are virtually the same. 

In much the same way we analyzed the homogeneous 
case in the previous section, here we obtain the equa- 
tion of state of the model along with the various ther- 
modynamic properties in the trapped system. We exam- 
ine the accuracy of the analytical approaches discussed 
in the previous section in this inhomogeneous setup and 
compare them against QMC simulations in a trap. (For 
the trapped bosons, however, the equations derived in 
the previous section cannot be reduced to simple an- 
alytical expressions as in the homogeneous case; they 



must be computed numerically.) In addition, we explore 
the regimes of validity of another approximate method, 
namely, the LDA which is based on the QMC results of 
the homogeneous system. 



urc ^h) shows the density in the center of the trap as 
a function of the characteristic density, a relation that is 
very useful to experimentalists for controlling their den- 
sity profiles by only changing the total filling in their sys- 
tems once the value of the trapping potential is known. 



A. Equation of state 

In a trapped system, where local densities change with 
position, the meaning of 'equation of state' may become 
somewhat unclear. This is because in an inhomogeneous 
trapped system the local phases are determined not only 
by the onsite interaction and total filling of the system 
but also through the strength of the confining potential. 
Fortunately, it turns out that the total filling and the 
curvature of the confining potential come into play only 
as a combination of the two: This is the characteristic 
density. It is the trap-equivalent of the overall density in 
homogeneous systems. 

The characteristic density (p) was introduced for lat- 
tice fermions in one dimension in Refs. [3l|, l32|, where 
it was shown by means of QMC simulations, that the 
scaled dimensionless variable p = pN[a/£)'^ can be used 
to generate a state diagram of the coexisting phases in 
the trap in the plane p vs U /t in a way that is inde- 
pendent of the specific values of the number of particles 
and the amplitude of the confining potential. Further 
studies validated this approach for the one-dimensional 
fermionic case within Bethe-Ansatz schemes J33l43a | , and 
state diagrams have been obtained for fermions in higher 
dimensions by means of dynamical mean-field theory 
(DMFT) [37[- The same idea has been shown to work 
with equal success in bosonic systems [l^ , in the frame- 
work of which the concept of characteristic density has 
been further discussed within the local density approxi- 
mation [sMi and finite-size scaling [3^, |4^ . It should be 
noted however that the concept of characteristic density 
remains a meaningful quantity even in cases where the 
LDA fails [13. 

As follows from the above discussion, one can then 
generate equations of state for trapped systems based 
on the characteristic density. This has recently been 
done by Roscilde in Ref. [30] for soft-core bosons in var- 
ious superlattice potentials. The equation of state for 
two-dimensional hardcore bosons is given in Fig. IHl^a), 
which shows the dependence of the characteristic density 
in our two-dimensional trapped systems in terms of the 
global chemical potential. The SSE results are denoted 
in the figure by circles, whereas the various approxima- 
tion schemes are represented by the different lines. The 
inset shows the deviations of the various approximation 
schemes from the QMC data. It is clear from the figure 
that, at least to some extent, all approximation schemes 
without exception provide a fairly accurate description 
of the equation of state. The spin-wave corrected results 
however provide a better match than the cluster-mean- 
field approximation which in turn shows only negligible 
improvement over the Gutzwillcr-mean-field results. Fig- 



24- 



a 13- 



-0.75 




/j/(2dt) 



FIG. 6: (Color online) (a) Equation of state and (b) density of 
bosons in the center of the trap vs p for harmonically trapped 
hardcore bosons in a two dimensional lattice (48 x 48 sites and 
V — 0.02). The circles indicate the QMC results whereas the 
dotted, dashed, solid, and dot-dashed lines are the Gutzwiller- 
mean-field, the cluster-mean-field, the Gutzwiller-mean-field 
with spin-wave corrections, and the LDA results, respectively. 
The inset in panel (a) shows the deviations of the various 
approximation schemes from the SSE data. 

Figure [S] also shows another useful approximation 
method that is usually applied to trapped systems. This 
is the LDA, which provides predictions of local quanti- 
ties based on results obtained from matching homoge- 
neous systems. Within the LDA, local observables of the 
confined system are approximated by their values in the 
corresponding homogeneous system, where the chemical 
potential of the homogeneous case is taken to be the cor- 
responding local chemical potential in the trap. In gen- 
eral, the LDA has been shown to give reasonably accurate 
descriptions of other confined systems, but is known to 
fail in regions where the local potential changes rapidly 
when compared with the correlation length of the rele- 
vant phase in the homogeneous case. The LDA results 
here were generated from the QMC results of the homo- 
geneous systems studied in the previous section. As the 
figure illustrates, the LDA provides a very good approx- 
imation of the QMC results for the density in the center 
of the trap, although, as expected, it underestimates the 
characteristic density required to form a Mott insulator 
in the center of the trap [Fig. ^h)] . 



B. Local densities and compressibilities 

In recent experiments, it has become possible to irnage 
local density profiles in optical lattices in-situ |4l| - |43| . 
This now allows experimentalists to explore the behavior 
of local quantities in a trap. Two quantities that provide 
insight into the properties of the trapped system are the 
local densities and compressibilities [l3|. In Fig. [3 we 
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show the local radial densities p - averaged over the az- 
imuthal angle - and the local compressibilities dp/d^Li 
both as a function of f^. The local compressibility was 
computed as a numerical derivative of the radial density 
when viewed as a function of the local chemical poten- 
tial, in much the same way as it would be derived in an 
actual experiment. This is done for two values of chem- 
ical potential: The top panels correspond to p ~ 4.62 
[p/{2dt) = 0.175] for which there is no Mott-insulating 
region in the center of the trap, while the bottom panels 
correspond to p « 15.02 [p/{2dt) = 1.3], in which case 
there is an insulating region in the center, as indicated 
by the plateau in the density profile when f is small. The 
QMC data presented in the figure correspond to two dif- 
ferent values of trap curvature which however trace the 
same curve, reflecting the fact that proper scaling elim- 
inates the dependence of the results on the strength of 
the trap. Interestingly, all approximation schemes are 
in good agreement with the exact QMC results. The 
spin- wave corrected results and the LDA provide the best 
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FIG. 7: (Color online) (a),(c) Density profiles and (b),(d) lo- 
cal compressibilities of harmonically trapped hardcore bosons 
in a two-dimensional lattice as a function of the normalized 
distance Vi from the center of the trap. The top panels cor- 
respond to a characteristic density p « 4.62 [p/{2dt) = 0.175] 
for which there is no Mott-insulating region in the center of 
the trap, while the bottom panels correspond to p ~ 15.02 
[p/{2dt) = 1.3], in which case there is an insulating re- 
gion in the center. The SSE results are indicated by cir- 
cles {V = 0.01) and squares {V — 0.02), whereas the dot- 
ted, dashed, solid, and dot-dashed lines correspond to the 
Gutzwiller-mean-field, the cluster-mean-field, the Gutzwiller- 
mean-field with spin-wave corrections, and the LDA results, 
respectively. The insets in the left panels show the deviations 
of the approximation schemes from the SSE results. 



match to the QMC results, although as is evident from 
the insets of the density profile panels, in the vicinity of 
the transition between the superfluid and the insulating 
phase, the LDA method shows larger deviations. This is 
also reflected in the compressibilities. 



C. Local phases 

The concept of characteristic density can also provide 
a general picture of the spatially-separated local phases 
inside the trap as the total filling and/or trapping po- 
tential are varied. To show that, we have computed the 
inner and outer radii fin and fout which enclose the su- 
perfluid region of the bosonic cloud. This was done by 
inspection of the local density profiles obtained via the 
QMC simulations for the different chemical potentials. 
In practice, the inner radius fin has been determined as 
the radius inside of which the local density is greater 
than 0.999 whereas fout was defined as the radius out- 
side of which the local density is lower than 0.001. We 
have verified that our results are robust against small 
deviations from the above values. Figure [5] shows the 
resulting diagram: Both radii are plotted for two differ- 
ent systems (two different trapping potentials) but with 
the same characteristic density at each point in the dia- 
gram. The circles indicate the SSE results with V = 0.01 
while squares correspond to y = 0.02. The lines in the 
figure indicate the predictions of the different approxi- 
mation schemes. This figure too illustrates that the con- 
cept of characteristic density is justified even when the 
LDA is not accurate, as the two systems corresponding to 
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FIG. 8: (Color online) Local phases of harmonically trapped 
hardcore bosons in a two-dimensional lattice; the notations 
E, SF and F correspond to empty, superfluid, and completely 
filled lattices, respectively. The SSE results are denoted by 
circles iy = 0.02) and squares {V = 0.01), whereas the 
various lines are the results of the different approximation 
schemes: Gutzwiller-mean-field (dotted lines), cluster-mean- 
field (dashed lines), Gutzwiller-mean-field with spin- wave cor- 
rections (solid lines), and LDA (dot-dashed lines). The sys- 
tems depicted here have 48 x 48 lattice sites. (Note that the 
Gutzwiller- and cluster-mean-field results here are concealed 
by the spin- wave-corrected results.) 
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the two values of trapping amplitudes trace virtually the 
same curve. Also, remarkably the various approximation 
schemes but the LDA yield very accurate predictions of 
both radii. 

Within the LDA, the two radii fjn and fout; at which 
the local densities arc {fii) = 1 and {fii) = 0, respectively, 
can be obtained in a straightforward manner: From the 
homogeneous results discussed in the previous sections, 
we know that these densities are reached at ij./{2dt) = 
±1. Translating back to the harmonic trap case, this 
condition becomes fi/{2dt) — rf /{2d) = ±1. Solving for 



for the two cases, wc arrive at: 



Vm/^ - 2d 



HnXDA = R.C 
r'out.LDA = \/ii/t + 2d. 



(37) 



The two radii are indicated by the dot-dashed lines plot- 
ted in Fig. m As one can see, the LDA yields a rather 
poor prediction of the two radii. This is expected how- 
ever since these radii mark the transition from the super- 
fluid phase to the insulating phase in the periodic system 
where diverging correlations are present. 
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FIG. 10: (Color online) Momentum distribution of harmon- 
ically trapped hardcore bosons in a two-dimensional lattice 
as a function of the radial momentum kr = \/k% -\- ky. (a) 
p f» 4.62 \p/{2dt) = 0.175] and (b) p « 15.02 [p/{2dt) = 1.3] 
for a system with 48 x 48 sites. The SSE results are indi- 
cated by circles whereas the lines are the results of the ap- 
proximation schemes: Gutzwiller mean-field (dotted lines), 
cluster-mean-field (dashed lines), and Gutzwiller mean-field 
with spin- wave corrections (solid lines). 



D. Thermodynamic properties and the momentum 
distribution function 



To conclude the analysis of the two dimensional case, 
here we present measurements of the kinetic and total 
(free) energies of the system. These are shown in Fig.[3]as 
a function of the characteristic density. The SSE results 
are denoted by circles, whereas the dotted, dashed and 
solid lines indicate the Gutzwiller-mean-field, the cluster- 
mcan-field, and the spin-wave-corrected approximations. 
Both the kinetic energy and the inset in the free-energy 
panel, indicate that spin-wave-corrected results provide 
a far better match than the two mean- field methods. 
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FIG. 9: (Color online) (a) Free energy and (b) kinetic energy 
of harmonically trapped hardcore bosons in a two-dimensional 
lattice as a function of the characteristic density. The circles 
indicate the SSE results (48 x 48 sites and V = 0.02), whereas 
the dotted, dashed, solid, and dot-dashed lines are the mean- 
field, the cluster-mean-field, the mean-field plus spin-waves, 
and the LDA results, respectively. The inset of the free-energy 
panel shows the deviation of the approximation schemes from 
the SSE results. 



Next, we analyze the momentum distribution of the 
model ri(k) as a function of the radial-momentum coordi- 
nate kr = \ k"^ + fcy (where, as with the radial densities, 

the azimuthal angle is averaged over) . This quantity can 
be directly probed in experiments with ultracold atomic 
gases via absorption imaging after timc-of- flight [iHsf. 
In the homogeneous case, due to the long-range decay of 
one-particle correlations in the compressible phase, the 
momentum distribution function has a delta peak sin- 
gularity at kr = 0. In general, in the insulating phase, 
the off-diagonal one-particle correlations decay exponen- 
tially (they are zero in our hardcore model), yielding 
a broad momentum distribution (completely flat in our 
case). When a harmonic trap is present, the compress- 
ible and incompressible insulating phases coexist, so the 
zero- momentum peak is a smoother function of kr . Fig- 
ure [10] shows the behavior of the momentum distribution 
for two values of the characteristic density: p k, 4.62 
[fi/{2dt) = 0.175], for which there is no insulating phase 
at the center of the trap, and p w 15.02 [p/(2dt) = 1.3] 
for which such a phase exists. Interestingly, it is the lat- 
ter value for which the zero-momentum peak is higher, 
despite the presence of the insulating phase in the center 
of the trap. This is because the system corresponding to 
the higher characteristic density has a higher number of 
bosons in the supcrfluid domain that surrounds the Mott 
insulator. 

Looking at the Gutzwiller- and cluster-mean-field ap- 
proximations, they both provide rather accurate predic- 
tions away from zero momentum but overestimate the 
peak by « 16%. The spin-wave corrected prediction is 
much better in both cases, giving only a « 5% error in 
the worst case. 
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V. TRAPPED HARDCORE BOSONS IN THREE 
DIMENSIONS 

In the previous section, we showed that the Gutzwiher- 
mean-field solution aheady provides a fairly good de- 
scription of the hardcore Bose-Hubbard model in the 
presence of a harmonic potential, where the deviations 
from the QMC data were shown to be rather small in 
most cases. This was true both for global observables 
such as the energies shown in Fig. [5J and also for local 
observables such as local densities and compressibilities 
shown in Fig. [T] When spin-wave corrections were taken 
into account, the results were improved to a point where 
most errors dropped to below 1%, yielding virtually exact 
results. 

The encouraging results of the previous sections sug- 
gest that in three dimensions the spin-wave corrected 
results should be an even better approximation of the 
exact solution. In light of this, in what follows we ex- 
plore the harmonically trapped three-dimensional hard- 
core Bose-Hubbard model using the Gutzwiller-mean- 
field approach with and without the addition of spin- 
wave corrections. The three-dimensional harmonically- 
confined system is a good example of a system that be- 
comes very demanding computationally if one wants to 
study the ground state for lattice sizes that are on the 
same order of magnitude as the ones realized experimen- 
tally. Having a large linear system size is critical when 
probing local observables that change with distance as is 
the case of the harmonic trap. 

Below we present the results of the three-dimensional 
analysis: The Gutzwiller- mean- field (dotted lines), the 
spin-wave-corrected (solid lines) and the LDA (dot- 
dashed lines) predictions of the equation of state (Fig. 
[TT|) . the boundaries of the local phases in the trap (Fig. 
[T^ and the energies (Fig. [T5)) of harmonically trapped 
hardcore bosons in three dimensions (for a 24 x 24 x 24 
lattice). These results can be used in current experiments 
to understand density profiles and other properties of the 
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FIG. 11: (Color online) (a) Equation of state and (b) density 
of bosons at the center of the trap for harmonically trapped 
hardcore bosons in a three-dimensional lattice (24 x 24 x 24 
sites and V = 0.08). The dotted lines indicate the Gutzwiller- 
mean-field while the solid and dot-dashed lines are the spin- 
wave-corrected and LDA results, respectively. 
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FIG. 12: (Color online) Local phases of harmonically trapped 
hardcore bosons in a three-dimensional lattice; the notations 
E, SF, and F correspond to empty, superfluid, and completely 
filled lattices, respectively. The Gutzwiller-mean-field results 
are denoted by the dotted line, whereas the dashed and solid 
lines correspond to LDA results and spin-wave corrected re- 
sults, respectively. (Note that the Gutzwiller-mean-field re- 
sults here are concealed by the spin-wave-corrected results.) 



system in regimes in which the local filling does not ex- 
ceed one and the ratio U/t is very large. 
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FIG. 13: (Color online) (a) Free energy and (b) kinetic en- 
ergy of harmonically trapped hardcore bosons in a three- 
dimensional lattice as a function of the characteristic density. 
The dotted lines indicate the Gutzwiller-mean-field results 
whereas the solid lines are the spin-wave-corrected results. 
Here, V = 0.08. 

As all three figures indicate, in three dimensions the 
spin-wave corrections are less pronounced than in two 
dimensions. However, they still provide a discernible im- 
provement over the mean-field results, particularly for 
the kinetic energy. 

In Fig. I14[ we examine some of the local properties 
of the three-dimensional system, specifically, the radial 
density profiles p{f) and the local compressibilities, as 
they are predicted by the Gutzwiller-mean-field, the spin- 
wave-corrected and the LDA methods. The figure shows 
the results for two values of chemical potential: The top 
panels correspond to p « 23.90 [fi/{2dt) = 0.3] for which 
there is no Mott-insulating region in the center of the 
trap, while the bottom panels correspond to p « 118.12 
[/i/(2dt) = 1.5], in which case there is an insulating region 
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in the center, as indicated by the plateau in the density 
profile when f is small. As is evident from the figures, 
in the vicinity of the transition between the superfluid 
phase and the surrounding insulating phases which are 
rather abrupt, the LDA results detach from the mean- 
field and spin- wave ones. This is also reflected in the 
compressibilities . 



2.5 



5.?:i 



0.6 



'^ 



■5 0.: 



„ 0.5 






0.05 




a. 



0.1 



0.05 



2.5 



5.0 



2.5 



FIG. 14: (Color online) (a),(c) Density profiles and (b),(d) lo- 
cal compressibilities of harmonically trapped hardcore bosons 
in a three-dimensional 24 x 24 x 24 lattice and {V = 0.08) 
as a function of the normalized distance fi from the center 
of the trap. The top panels correspond to a characteristic 
density p ^ 23.90 [ii/{2dt) = 0.3] for which there is no Mott- 
insulating region in the center of the trap, while the bottom 
panels correspond to p ~ 118.12 [^/{2dt) = 1.5], in which case 
there is an insulating region in the center. The Gutzwiller- 
mean-field results are indicated by the dotted lines, whereas 
the dashed, solid, and dot-dashed lines correspond to the 
Gutzwiller- mean-field, the cluster-mean-field, the Gutzwiller- 
mean-field with spin-wave corrections, and the LDA results, 
respectively. 



VI. SUMMARY AND CONCLUSIONS 

In this paper, we have studied the ground-state proper- 
ties of homogeneous and trapped inhomogeneous systems 
of hardcore bosons in two- and three-dimensional lat- 
tices. We have obtained results both by means of quan- 
tum Monte Carlo (SSE) simulations, and via several an- 
alytical approximation schemes, such as the Gutzwillcr- 
mean-field approach, a cluster-mean-ficld method, and a 



spin- wave analysis which takes quantum fluctuations into 
account. We studied the equation of state and various 
quantities of interest such as densities and compressibil- 
ities, condensate and superfluid densities, as well as the 
momentum distribution function and other basic thermo- 
dynamic properties. 

In the homogeneous case, we showed that the equation 
of state can be predicted reasonably-well by all approx- 
imation schemes. However, for other quantities such as 
the kinetic energy, and the superfluid and condensate 
densities, the Gutzwiller and cluster mean-field predic- 
tions were found to deviate considerably from the quan- 
tum Monte Carlo results, particularly close to half fill- 
ing. Interestingly, the addition of quantum fluctuations 
through spin-wave corrections was shown to yield virtu- 
ally exact results for all quantities but the compressibil- 
ity. 

In the inhomogeneous case, we found that for the equa- 
tion of state and some quantities, such as the onsite 
densities, even the simplest method, i.e., the Gutzwiller- 
nrcan-field approach, is capable of providing a fair an- 
alytical description. For other quantities, on the other 
hand, Gutzwiller-mean-field theory is not a very good 
approximation. For example, for the momentum distri- 
bution function we found that it largely overestimates 
the very low momenta occupations. A cluster-mean-field 
approach, based on larger cell sizes, did not prove to be 
much of an improvement over the Gutzwiller-mean-field, 
especially considering the large amount of additional de- 
grees of freedom in the ansatz. 

The addition of spin-wave corrections turned out to 
be an excellent approximation for all quantities, yielding 
virtually exact results when compared against the QMC 
data. This, in turn, allowed us to study large systems of 
three-dimensional harmonically-trapped hardcore bosons 
of sizes that would be very demanding for QMC tech- 
niques to simulate, with a fraction of the effort. Hence, 
this analytical approach can be used to study trapped 
systems in which the onsite repulsion is much larger than 
the bandwidth and the filling-per-site that is always lower 
or equal to one. 

Finally, we examined the validity of the local den- 
sity approximation in these harmonically trapped sys- 
tems. We showed that it provides a fairly good de- 
scription of local quantities as long as the measurement 
takes place away from the superfluid-insulator bound- 
aries, which correspond to a phase transition in the ho- 
mogeneous case. 
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